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INTRODUCTION 

In the last decade, with the advent of high speed com- 
puters, numerical calculations using finite-differencing 
techniques applied to partial differential equations has 
been of great interest to gas dynamicists and meteorologists. 
This has been brought about by the accuracy, capacity and 
speed that electronic computers offer in the solution of the 
complex partial differential equations describing the motion 
of fluids. 

Finite-difference equations may be constructed and used 
in various ways depending on accuracy, stability and import- 
ant physical considerations, e.g., conservation laws. A 
"conservation law form" of a system of differential equations 
may take the form 

E + 1^- F = ° (1) 

where F and E are conservative variables. Another form of 
equation (1) may also be considered. This is the "advective 
form" of equation (1) which is 

5t f: + M h Fl = 0 (2) 

where F, and E are vectors and R is a matrix. 
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In development and use of finite-differencing techniques, 
meteorologists have favored the advective form of the partial 
differential equations. This is partly due to the form of 
the advection equations which are encountered in meteorological 
studies. Molenkamp (9) and Crowley (4) have extensively ap- 
plied advective differencing methods to meteorological model 
problems and have examined the results. However, they have 
not been able to obtain accurate solutions to the problems 
considered, except when relatively complicated and time-con- 
suming fourth-order accurate techniques were employed. On 
the other hand, gas dynamicists have particularly been 
interested in the application of differencing techniques to 
"conservative law form" of partial differential equations. 
Kutler (5) and Anderson and Vogel (1) have successfully 
applied conservative techniques to sonic-edged, conical, 
wing-body combinations and flow about a rectangular wing 
moving supersonically. 

It is the purpose of this research to investigate the 
possibility of application of the most recent conservative 
and widely applied numerical techniques in gas dynamics to 
problems encountered in meteorological computations. 

Solutions of the advection equation are obtained using con- 
servative differencing methods common to gas dynamics. 

These results are compared to those obtained by Molenkamp, 
who differenced the advection equation directly. The 
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comparison shows that better results are obtained where 
conservative form of governing equation is used. This is 
in agreement with the results obtained by Crowley (4) . In 
addition, better results are obtained with lower order con- 
servative methods as compared with higher order differencing 
applied to the advection equation. 
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PROBLEM DEFINITION 

In a two-dimensional rectangular coordinate system, the 
differential advection equation is 


3 A 

3T 


3A 

37 


3 A 


(3) 


where A is the quantity being advected. Velocity components 
u and w are respectively in the x and z directions and t is 
time. If a steady velocity field is chosen then velocity 
components are no longer functions of time. The above 
equation becomes linear and an analytical solution is then 
possible. Considering the motion of the fluid to be a rotation 
with constant angular velocity, ft, an equivalent form of 
equation (3) in cylindrical coordinates becomes 


3 A 
Tt 


+ ft 


3A 


0 


(4) 


when 0 is the angular coordinate and radial velocity is zero. 
Equation (4) is the wave equation and its analytical solution 
is found to be 

A(r,9,t) = Ao(r,0-ftt) (5) 

where r is the radial distance from the axis of rotation and 
Ao is the initially given distribution of A at time zero. 
Equation (5) shows that the solution of the wave equation (4) 
is an angular displacement of the initial distribution Ao . 
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A conservative form of equation (3) may be obtained con- 
sidering an incompressible flow. Then the continuity equa- 
tion is 


or 


V.q = 0 


3u 3w 

ax az 


= o 


(6) 


(7) 


multiplying each side by A 


A + A If “ 0 


(8) 


Now add equation (8) to equation (3) , or 


aA . a a 3 a 

3t + u 53r + w 3? + 


au 

31 ? 


3w 

37 


= 0 


(9) 


which may be written in the following form 


TT7 + 37 (Au) + (Aw) = 0 


( 10 ) 


This is the conservative form of equation (3). The differ- 
encing techniques that are described in the next section are 
applied to the general form of this equation. 
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NUMERICAL TECHNIQUES 


In the following section those differencing methods con- 
sidered in the present paper are explained. The form of each 
equation is given when applied to the general conservative 
hyperbolic partial differential equation in two dimensions. 


3E DF 3G 
7t + -3x + 3? ~ 


0 


( 11 ) 


The modified equation for each technique is obtained by apply- 
ing these methods to the one-dimensional wave equation 




0 


( 12 ) 


In addition to the above, the stability criterion for each 
technique is given as obtained from the linear stability 
analysis. 


Brailovskaya Method 

The first-order predictor-corrector scheme described 
below was devised by I. Y. Brailovskaya (2) based on central 
differencing. When Brailovskaya's technique is applied to 
equation (11) , the result is 


=-n+l p n At , n n 

j,k j,k 1 Ax ' j+l,k j-lfk 


At , n _ n . 
THy G j ,k+l G j,k-1 } 
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„n+l _ „n _ At ,«n+l _ «n+l > 

E j,k ~ j,k 2 Ax j+l,k j-l,k 


At /*-n+l jsrn+l . 

TKy tG j,k+l ‘ G j ,k-l 


(13) 


The modified partial differential equation (6) for 
Brailovskaya's method may be found by applying this scheme 
to one-dimensional wave equation (12) 


u, + cu = w cv Ax u 


XX 


( 14 ) 


where v is the Courant number. 

Brailovskaya's technique is easy to program because of 
the simplicity of the structure of the scheme and the simil- 
arity of the differences in both predictor and corrector. 

The latter allows the programmer to define only one set of 
boundary conditions for both of the above steps. 

This first-order technique is stable under the follow- 
ing conditions. 


o At 


1 C At 

max Ax 

$ 1 , 

| max Ay 


where cr max is the maximum eigenvalue of the hyperbolic sys- 
tem under consideration. 

For better accuracy of the computation, the modified 
equation above requires the mesh size to be small such that 
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it decreases the magnitude of the second-order error. 

Brailovskaya's technique is not widely used because of 
its low order of accuracy and also because of the predictor- 
corrector sequence form which increases the computation time 
to that of the second order techniques. 

Lax-Wendroff Method 

A second-order differencing scheme was derived by Lax 
and Wendroff (7) for which the stability criterion is defined 
by 


At 


At 

0 -r— 

max Ax 

^ 1 / 

u max Ay 


where again is the maximum eigenvalue of the hyperbolic 

system under consideration. This technique when applied to 
equation (11) yields 


n+1 n _ At , p n n . At . n _ n 

j,k j,k 2Ax j+l/k j-l,k 2Ay ' j,k+l j,k-l 


+ k (tkk.) 

2 'Ax' 


A'" (F n . -F n )-A' n - (F n , ■ -F n 1 ) 

j4,k ^ + 1 ' k 3' k j-±,k =>' k ^ 1 ' k 


+ k (—) 

2 Ay ^ 


B' n . (G n -G n ,)-B' n (G n ,-G n ) 
j,k+J l^tl 3 ,k jk .l 3,k D,k-1 


, 1 (At) 


8 Ax. Ay 


a " n ic n _c n i-A' n -r n 

j + lA j+l,k+l j+l,k-l' j-l,k+l j-l, k 


3 ,k+l u j+l, k+1 j-1 ,k+l' j,k-l v 3+l,k-l j-l,k-l 


( 17 ) 
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where 


and 




A"" (E) 

j+ 7 , k 


A' 



n 

j + lA 


+E 


n 

jr* 


) 


(18) 


B" n (E) 

j+?,k 


B ' [7 {E j+l / k +E j f k ) .. 


(19) 


The Lax-Wendroff method applied to the one-dimensional wave 
equation forms the following modified partial differential 
equation 


u fc + cu x = - ^ c (1 -v 2 )Ax 2 u xxx + (20) 

The form of the above equation confirms the order of accuracy 
of the Lax-Wendroff technique. One may note that at a Courant 
number equal to unity, the above equation reduces to the exact 
wave equation (12) and thus provides an exact solution. 

MacCormack Method 

MacCormack (8) developed a second-order predictor- 
corrector sequence for use in studies involving hypervelocity 
impact cratering. When applied to equation (11), it yields 
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TrfW-1 - V n _ At fv n -v n At ir n -r n ) 

E j,k “ E j,k 2x ^ F j+l,k F j,k ffy j,k+l j,k 


= 1 [e!? ' +g n+ i-^(F n+ f-Ff^ . 

j,k 2 j /k j,k Ax 3 ,k j-1,1 


At ,*n+l »n+l . 

- Zy ( 5 j,k- B j,k-l>J 


It is interesting to note that this technique is a prefer- 
ential scheme using a forward predictor and backward cor- 
rector. The backward predictor and forward corrector version 
of MacCormack ' s technique is also examined and results are 
reported in this paper. MacCormack' s differencing scheme 
has been applied to gas dynamic problems in recent years 
and has resulted in accurate solutions comparable to better 
second-order methods (1, 5) . 

In this case, the stability bound is again found to be 


a 7T- < 1 
max Ax 


o TT— < 1 

max Ay 


The following modified partial differential equation is 
obtained when MacCormack' s technique is applied to equation 
(12 ) 



Note again that at Courant number of unity the above equation 
reduces to equation (12) . 
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Rusanov Method 

In 1969, Rusanov (10) and Burstein and Mirin (3) 
separately developed a third-order accurate scheme which 
has been of great interest to gas dynamicists where high- 
speed computers are available. This three-level predictor- 
corrector technique provides accurate solutions when applied 
to gas dynamic equations. The Rusanov method applied to the 
general equation (11) results in 


E 


(D = 1 

3 +|' k 4 7 


E "+i,k+i +E 5,K+ 1 +E 5+i,k +E 3,ic 


1 /At 
6 ^Ax 


j+1 , k+1 j + l,k j,k+l j , k 


At 

Ay 


4-P^ 

U j + l,k+l +Cj j,k+l ^j + l,k J,k 


(2) _ n 1 f A t 
j,k 3 ,k " 3 l Ax 


F d) +F U> - F d) - F U> 

.j.1 u.i .^1,1.1 ,1 .1,1 

j+j,k+j 3-jrk+j 



r (l) +r (l) r (l) r U) 

•4.1 -i i ^r G .^i k r. u i 

3+j'k+j 3-jfk+j J+jrk-j 3-jrk-j 


} 


n+1 ,n 1 f At ,„n _ p n At , p n _,„n 

j,k j,k ” 4 7 ^Ax F j+1 ,k j-l,k 37Tx j+2,k j+l,k 


+ 2F 5-l,k- F ?-2,k» + 5§< G ”,k+l- G ?,k-l> 


r~(C n , ,--2G n , ,.+2G n , .-G n . J } 
>Ay J , k+2 3, k+1 3,k-l 3,k-2' 
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_ 3 .At (2) _p (2 ) At (2) (2) n 

8 l A^ (F j+l,k F j-l,k )+ Sy (G j,k+l G j ,k-l 

+ ^30 j+2 , k” 4E j+1 /k +6E j , k -4E j-1 , k +E j-2 , k } 

+ Y 30 ^; /k+ 2- 4E ?,ktl^,k- 4E ;,k-l^,k-2 } 


This technique is stable when 

At 


°max Ax 


$ 1 


a 


max 


At 

Ay 


< 1 


and 


where 


2 4 

4v -v $ - a) $ 3 . 0 


w = - 24 y 


30 


When Rusanov's technique is applied to equation (12), the 
following modified partial differential equation results. 


u t + cu x = 


1 A 3 

'27 c Ax 


u> . . 3 

— -4v + v u 
v xxxx 


1 . 4 

m c Ax 


5u-4-15v^+4v 4 ] u 


xxxxx 


Again it may be noted that when 


v = 1 . 0 


w = 3.0 


an exact solution of wave equation is formed. 


(24) 

(25) 

(26) 

(27) 

(28) 
(29) 
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Kutler-Warming Method 

The most recent third-order differencing scheme developed 
is the modified version of the Rusanov three-level predictor- 
corrector sequence developed by Kutler and Warming (11) . 

When applied to equation (11) , the Kutler-Warming technique 
yields 


E (1) ~ 

E< 2 > = 
3/k 

n+l 

E j.k 


n 2 At ,_n _n ,2 At, n _ n . 

j,k 3 Ax j + l,k j,k T Ay j,k+l j,k 


E n At (F U)_ P (l> , !lt uuai 

3,k ],k 3 Ax ],k ]-l,k 3 Ay 3,k 3,k-l 


E n + Y 

E j,k + Y 30 


41 — 4E 

j+2 ,k 4E 

+ Y-^nfE 1 ? -4P n 

Y 30[ £ 'j, 


,n +6E n — 4E° +E n 
'j+l,k +6E j,k 4E j-l,k +E j-2, 

,. 0 -4E n ,_+6E n . -4E n . .+E n . 
k+2 3,k+l 3,k 3,k-l T,k- 

- -mjr f-2F n _ , +7F n -7F n .+2F 1 ? . 

24 Ax L 3 + 2, k 3 + 1 / k 3 ~l,k 3-2, k, 

- f-2G n +7G n lr . 1 -7G n . n +2G n . , 

24Ay [ 3 , k+2 3 ,k+l 3 ,k-l 3,k-2 

3 At f (2) (2) 

8 Ax L j+l/K 3-1 


' _3 At 

r r (2) _ c (2) 

8 Ay 

[ G j ,k+l G j ,k-l. 


The modified equation which results from the application of 
this technique to the linear wave equation (12) is 


u t +cu x = - IT c iJ<3 


[S -4u + v 3 ] u 


XXX 


1 A 4 

m c4x 


5w-4-15v 2 +4v 4 


u 


XXXXX 


(31) 
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This technique is modified so that results comparable 
with the Rusanov scheme are obtained in less computation 
time. Similarity of the results come from the fact that 
their modified equations are similar, and, in fact, identical 
up to the fourth-order. A decrease in computation time is 
the direct result of simpler equations and less computation 
in all three levels due to elimination of intermediate grid 
point calculations. Note again when 


v = 1.0 


to = 3.0 


(32) 


the modified equation reduces to equation (12). The usual 
stability requirement is 


, At 

max Ax 


<1 


max Ay 


$ 1 


(33) 


and 


2 4 

4v - v < to < 3.0 


( 34 ) 
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SAMPLE PROBLEM 

An initially specified disturbance is considered in 
rotation about an axis normal to the XZ-plane or the plane 
of grid points. The angular velocity ft is taken to be 
constant for which the stream function is defined as 



(35) 


where 


r = 


(x-x") + 


(2-2") 2 


(36) 


and x' and z' are the coordinates of the intersection point 
of the axis of rotation and the grid plane. 

An incompressible flow is considered where the govern- 
ing equation of the fluid flow is 

3A , 3 (Au) , 3 (Aw) _ „ 

3t 3x 3z " U 1 


Velocity components u and w are determined from the follow- 
ing equations 




(38) 


The initial distribution of A is chosen as 


A 0 (x, z) 



for r"< 4A 


0 


for r"^ 4 A 


(39) 
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where A is the grid interval and 

(40) 

where x M and z" are the coordinates of the point of the maxi- 
mum value of A. 

It may be noted that the prescribed initial distribution 
above describes a cone with its base on the XZ-plane. The 
peak is placed such that the disturbance is away from the 
axis of rotation and does not hit the boundaries of the grid 
plane in the course of its rotation about (x' ,z') . Also, 
the values of A are forced to be equal to zero along the 
boundaries during the sequence of numerical integration. 

The following constants were used in all the calcula- 
tions except where otherwise specified. 

n = number of time iterations = 40 

= angular velocity = -0.001 radian/second 
At = time interval = 30 seconds 
Ax = Az = A = 1.0 
x ' = 12A 
z' = 12A 
x" = 18A 
z" = 12A 

All of the computations were performed on an IBM 360-65 
digital computer at the Iowa State University Computation 


" = [_ (x— x " 


il 


> + (z-z") 


Center . 
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RESULTS AND CONCLUSION 

The numerical techniques described earlier were first 
applied to the modified Burger's equation in one-dimension 

Conservative : 

Advective: 

An initial distribution of 


u (x ) = 1.0 

for 

0 

$ x $ 

50 

u (x) =0.0 

for 

50 

< X $ 

100 


was assumed. The results obtained were in complete agree- 
ment with Crowley's results indicating that conservative 
techniques are to be preferred over advective methods. 

Another important conclusion is also derived. Conservative 
methods are preferred for problems with continuous and smooth 
solutions as shown by Crowley and also in problems involving 
discontinuities such as shock waves. In fact, non-conservative 
differencing can result in improper wave speed in flows in- 
volving discontinuities. 

In order to further examine the conservative differenc- 
ing techniques described earlier, they were applied to the 
sample problem discussed in the last section with governing 
equation (37) . The relationship of the solution obtained by 
difference approximations to the analytic solution may be 


+ JL< U \ - n 

Jt + 3x 7 " 0 


3u 

Tt 


3u 

*3x 
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better represented when contour plots of A are examined. 

Also, the accuracy of each solution in terms of its general 
approximation, phase and radial displacement can be easily 
investigated (Figures 1-15) . In these figures, the analytical 
solutions are represented by broken lines, while the solid 
lines are the solutions obtained by differencing techniques. 
Table 1 contains a summary of several important features 
of the results obtained by Molenkamp, while Table 2 contains 
the results of present conservative solutions. 

Examining Figures 1, 2 and 3 with results given in 
Tables 1 and 2 for the first-order schemes, the superiority 
of conservative approximations is clearly shown over advective 
solutions. The advantage is more obvious when higher values 
of A isolines are considered. Computation time is 30% higher 
for Brailovskaya's technique than either Upstream N or 
Upstream N+l, but its higher maximum isoline approximation 
and lower radial displacement justifies its use. It is to be 
noted that the accuracy of these first-order techniques may 
be increased somewhat by decreasing the mesh ratio, — and 

AX 

. However, first-order techniques are only simple means 
of determination of the general behaviour of the solutions 
and therefore are not recommended for use in solution of 
complicated partial differential equations. 

Second-order techniques (Figures 4-9) resulted in 
generally better solutions than the first-order methods as 
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expected. The quality of the approximations of A isolines, 
using MacCormack and Lax-Wendroff conservative schemes are 
comparable to those approximations obtained by Molenkamp 
using advective Leap-Frog, Arakawa-Euler , and Arakawa-Adams- 
Bashforth techniques with conservative angular displacement 
error being 16 to 66 percent less than errors involved in 
advective solutions. 

In general, the MacCormack differencing scheme is a 
better method overall than any other second-order advective 
or conservative technique considering the general approxima- 
tion, computation time, error, and structure of the differ- 
encing equations. 

Third-order techniques (Figures 10-15) resulted in the 
most accurate solutions obtained in this investigation. 

The accuracy of Rusanov-Burstein-Mirin technique had been 
investigated (1) where accurate solutions were obtained for 
gas dynamic model equation, i.e.. Burger's Equation. The 
Kutler-Warming method is basically a modified form of the 
Rusanov technique and, in fact, the similarity of their 
modified equations suggests a close agreement of the solutions 
This proved to be true for the problem under consideration in 
this report. The approximation obtained by application of 
these third-order methods closely follows the circular pat- 
tern of the analytical solutions proving their advantage over . 
any lower order advective or conservative scheme. Tests 
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were made for different combinations of values of y^q and 
At. When the values 

y 30 = " Sff ' At = 60 

are used, the Kutler-Warming technique yields the best 
solution to the above problem. It is interesting to note 
that the above values of y^ and At correspond to about 
one-seventh of the lower bound for y^Q prescribed by equation 
(34) . Anderson and Vogel also found that better results are 
obtained when the y^g values corresponding to lower bound and 
fractions of the lower bound of the stability equation (34) 
were used (11) . This indicates that linear stability analysis 
resulting in equation (34) does not define accurate stability 
bounds for all linear and non-linear problems. It should be 
noted that y^ Q may not assume the value of zero, and there- 
fore a limit exists on how small the value of j y | is to be 
chosen. This is also shown in Figures 10-15 and Table 2. 

The quality of the solution is degraded as y 3Q assumes 
values lower than ~ . In general, the Kutler-Warming tech- 
nique is preferred over the Rusanov-Burstein-Mirin method 
mainly because of the simpler structure of the differencing 
equations in all three levels and elimination of intermediate 
grid calculations. A direct result of this is a considerable 
decrease in computation time. A comparison of the approxima- 
tions obtained from the application of the above conservative 
third-order schemes with the results of the advective 
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Roberts-Weiss method, indicates again, the superiority of 
the conservative differencing over advective computation. 
Comparable approximations were obtained in both cases with a 
ratio of 1.2 to 45 of required computation time in favor of 
conservative differencing (Tables 1 and 2) . 

In general, conservative differencing is to be pre- 
ferred over the advective approximation. Numerical experi- 
ments by Crowley and those reported in this paper confirm this 
fact. An improvement in results is obtained when higher order 
differencing techniques are applied. This is shown to be 
independent of whether the equation is in conservative or 
advective form. Those differencing techniques discussed in 
this paper are mainly gas dynamic differencing methods, 
but the results of this investigation in comparison with those 
by Molenkamp and Crowley indicate that these techniques may 
be applied to advection equations as well, resulting in 
better accuracy and more economical computation. The Kutler- 
Warming version of Rusanov's third-order technique resulted 
in the most accurate solutions and along with its short 
computation time presents, at the present time, an optimum 
differencing method in the solution of meteorological and 
gas dynamics equations. 



Table 1. Maxima, lag, and computation time after 40 iterations (advective) 
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Figure 6. 


Solution using the MacCormack (forward predictor, 
backward corrector) Method. 





Figure 7. Solution using the MacCormack (forward predictor, 
backward corrector) Method. 



Figure 8. Solution using the MacCormack (backward predictor, 
forward corrector) Method. 





Figure 9. Solution using the MacCormack (backward predictor, 
forward corrector) Method. 
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Figure 14. Solution using the Kutler-Warming Method. 
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